Hydration of Methanol in Water 
A DFT-based Molecular Dynamics Study 
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We studied the hydration of a single methanol molecule 
in aqueous solution by first-principle DFT-based molecular 
dynamics simulation. The calculations show that the local 
structural and short-time dynamical properties of the wa- 
ter molecules remain almost unchanged by the presence of 
the methanol, confirming the observation from recent exper- 
imental structural data for dilute solutions. We also see, in 
accordance with this experimental work, a distinct shell of 
water molecules that consists of about 15 molecules. We 
found no evidence for a strong tangential ordering of the water 
molecules in the first hydration shell. 

INTRODUCTION 

The solvation of alcohols in water has been studied 
extensively, [jjj It is of fundamental interest in physics, 
chemistry and biology, but also of importance in technical 
applications. The characteristic hydroxyl group allows 
alcohols to form hydrogen bonds and is responsible for 
the good solubility of the smaller alcohols. In contrast, 
the alkyl group is hydrophobic and does not participate 
in the hydrogen bonding network of water. The presence 
of both hydrophobic and hydrophilic groups make the 
microscopic picture of solvation of alcohol in water a non- 
trivial and therefore interesting matter. 

Understanding the solvation of methanol in water is 
a prerequisite for the study of chemistry of alcohols in 
aqueous solution. Important examples of such reactions 
are the conversion of ethanol into acetaldehyde in biolog- 
ical systems or the industrial ethanol production by acid- 
catalysed hydration of ethylene. An accurate microscopic 
understanding of the mechanism and kinetics of such re- 
actions is of fundamental interest. However, presently, 
this picture is still far from complete. Density Functional 
Theory (DFT) based Molecular Dynamics simulation has 
proved to be a promising tool provide such an insight. An 
accurate calculation of the chemical bonding is incorpo- 
rated via a DFT-based electronic structure calculations. 
The effect of temperature and solvent on the reactive 
events is implicitly accounted for via the Molecular Dy- 
namics technique. The implementation of DFT-based 
MD as proposed by Car and Parrinello j^j has proven to 
be extremely efficient. It has successfully been applied 
to study of a large variety of condensed-phase systems 
at finite temperature. Applications to chemical reactions 
include the cat-ionic polymerization of 1,2,5-trioxane ||, 
or the acid-catalysed hydration of formaldehyde Q . 

As a first step towards the study of chemical reac- 
tions involving alcohols we present in this paper a Car- 
Parrinello Molecular Dynamics (CPMD) study of the hy- 
dration of the simplest alcohol (methanol) in aqueous 



solution. Recent experimental work || has provided de- 
tailed structural information on the solvation shell. Var- 
ious molecular simulation studies (e.g. Ref. J6| — |lOt] have 
addressed structure and dynamics of both the solute and 
the solvent. This experimental and numerical work has 
revealed that there is a distinct solvation shell around 
the methanol, and that the water structure is little af- 
fected by the presence of a methanol molecule. In this 
paper we will address these structural properties and in 
addition consider the dynamics of the methanol and the 
water molecules in the solvation shell. 

This paper is organized as follows. First we outline 
the computational approach and its validation. Then we 
present the results for the structure and dynamics of a 
single solvated methanol in water. We conclude the paper 
with a summary and discussion. 

METHODS AND VALIDATION 

Electronic structure calculations are performed using 
the Kohn-Sham formulation of DFT. Q We em- 
ployed the BLYP functional, that combines a gradient- 
corrected term for the correlation energy as proposed by 



Lee, Yang and Parr |14| with the gradient correction for 
the exchange energy due to Becke |l3| . Among the avail- 
able functionals, the BLYP functional has proven to give 
the best description of the structure and dynamics of wa- 
ter. 15 1(| All calculations Jl7j were performed using the 
CPMD package, ff 

The pseudopotential method is used to restrict the 
number of electronic states to those of the valence elec- 
trons. The interaction with the core electrons is taken 
into account using semi-local norm-conserving Martins- 
Troullier pseudopotentials. Jl!| The pseudopotential cut- 
off radius for the H was chosen 0.50 au. For O and C 
the radii are taken 1.11 and 1.23 a.u. for both the l=s 
and l=p term. The Kohn-Sham states are expanded in a 
plane-wave basis set matching the periodicity of the peri- 
odic box with waves up to a kinetic energy of 70 Ry. Test 
calculations showed that for this structural and energetic 
properties were converged within 0.01 A and 1 kJ/mol, 
respectively. Frequencies are converged within 1 %, ex- 
pect for CO and OH stretch modes that are underesti- 
mated by 3 % and 5 % compared to basis-set limit values. 

To validate the computational methods outlined above 
we performed a series of reference calculations of rel- 
evant gas-phase compounds with the CPMD package. 
Energetics and geometry were calculated for methanol, 
water, two mono-hydrate configurations, and the di- 
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hydrate configuration shown in Fig. These calcula- 
tions were performed using a a large periodic box of 
size 10x10x10 A 3 . The interactions among the periodic 
images were eliminated by a screening technique sim- 
ilar to that of Ref. |20|| . In addition we determined 
for the methanol molecule both the harmonic vibra- 
tional frequencies and the frequencies at finite temper- 
ature (T= 200 K). The latter includes the anharmonic 
contributions, and were obtained from the spectrum of 
the velocity auto correlation function (VACF) of a 3 ps 
CPMD calculation at E= 200 K. The calculated peak po- 
sitions can be compared with experimental spectra. Re- 
sults of the gas-phase calculations were compared with 
results obtained with a state-of-the-art atomic-orbital 
based DFT package (ADF and with results from 

MP2 calculations of Ref. p2 1. In the comparison of the 
energies zero-point energies were not taken into account. 

(a) 




FIG. 1. Energy-optimized geometries of 

two water/methanol dimers and a trimer. Distances (A) and 
angles (degrees) are shown for three computational methods: 
CPMD-BLYP (top, present work), ADF-BLYP §l| (middle, 
present work) and MP2 Q (bottom). 

Complexation energies and geometries of the methanol 
hydrates are given in Tab. | and Fig. |l|. Deviations among 
CPMD and ADF are within 1 kcal/mole for the energies, 
smaller than 0.005 A for the inter-molecular bonds and 
within 0.03 A for the weaker intra-molecular bonds. This 
indicates a state-of-the art accuracy for electronic struc- 
ture methods employed in CPMD. Differences among 
BLYP and MP2 are within acceptable limits, with BLYP 
complexation energies smaller by 4 kJ /mole (dimer) and 
10 kJ/mole (trimer). These deviations are similar to the 
comparison of BLYP and MP2 for the water dimer bind- 
ing energy, (lj]]2j| where BLYP is 4 kJ/mole smaller, 
with the MP2 energy only 1 kJ/mol below the experi- 
mental value. Assuming similar differences for the com- 



plexation energies bonds in the methanol hydrates would 
suggest that BLYP underestimates the methanol-water 
binding energy by approximately 5 kJ/mol. Inter- and 
intra-molecular BLYP bond lengths are up to 0.02 and 
0.06 A longer compared to the MP 2 results, respectively. 

TABLE I. Complexation energies (kJ/mol) of methanol 
hydrates shown in Fig. jjj. Numbers are bare values without 
zero-point energy corrections and entropy contributions. 



Complex 


CPMD-BLYP 


ADF-BLYP° 


MP2 6 


CH 3 + H 2 (a) 


20.2 


20.2 


24.4 


CH3O + H a O (b) 


17.1 


17.6 


21.0 


CH3O + 2 H 2 


58.3 


59.6 


68.8 



a Ref. |32|. 

6 G2(MP2) method. MP2(full)/6-311+G(d,p) optimized ge- 
ometries. From Ref. [E2 1. 



Vibrational frequencies are listed in Tab. || A gam 
comparison of CPMD and ADF is excellent, consis- 
tent with the results for the energetics and geometries. 
Comparing the calculated finite-temperature frequencies 
against the experimental values shows that BLYP tends 
to underestimate the frequencies of almost all modes by 
« 10 %. This trend is a known feature of BLYP. For 
example similar deviations are observed for BLYP calcu- 
lation of water. [|l5| 

Overall we conclude that the reference calculations of 
gas-phase provides confidence that DFT-BLYP performs 
with a sufficient accuracy for a quantitative study of 
methanol hydration. 



TABLE II. Harmonic and T=200 K vibrational frequen- 
cies of gas-phase methanol molecule. 



mode 


Harmonic 
v (cm -1 ) 
CPMD-BLYP ADF-BLYP" 


Anharmonic 
v (cm -1 ) 
CPMD-BLYP Exp. 6 
(T=200 K) 


r(OH) 


280 


380 


280 


270 


v\cO) 


940 


950 


880 


1034 


r(CH 3 ) 


1040 


1050 


980 


1075 


r(CH 3 ) 


1130 


1130 


1070 


1145 


<5(OH) 


1330 


1340 


1270 


1340 


<5(CH 3 ) 


1430 


1430 


1320-1430 c 


1454 


<5(CH 3 ) 


1460 


1460 


1320-1430 c 


1465 


5(CH 3 ) 


1470 


1470 


1320-1430 c 


1480 


«/(CH 3 ) 


2940 


2910 


2640 


2844 


f(CH 3 ) 


2990 


2950 


2740 


2970 


f(CH 3 ) 


3060 


3020 


2830 


2999 


//(OH) 


3550 


3590 


3310 


3682 



Ref. 1 321. 
Ref. |§. 

Modes not separated. 



Broad peak with width listed. 
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SOLVATION 

We performed Car-Parrinello Molecular Dynamics 
simulations of the solvation of a single methanol 
molecule. We considered two systems: one with 31 water 
molecules and the other with 63 water molecules, yield- 
ing methanol- water solutions with mole ratios of 1:31 
and 1:63. In the following they are referred to as the 
small and large system, respectively. For reference we 
also performed a simulation of a pure water sample of 
32 molecules. The molecules are placed in a periodic cu- 
bic box with edges of 9.98 A (small solvated methanol 
system), 12.50 A (large solvated methanol system), and 
9.86 A (pure water) corresponding to the experimental 
densities at ambient conditions. The temperature of the 
ions is fixed at 300 K using a Nose-Hoover thermostat 
pi] p6| . The fictitious mass associated with the plane- 
wave coefficients is chosen at 900 a.u., which allowed for 
a time step in the numerical integration of the equations- 
of-motion of 0.145 fs. The two systems were equilibrated 
for 1 ps from an initial configuration obtained by a force- 
field simulation. Subsequently we gathered statistical av- 
erages from a 10 ps trajectory of the 31+1 molecule sys- 
tem, from a 7 ps trajectory of the 63+1 molecule system, 
and from a 10 ps trajectory of the pure water system. 

Structure 

In Fig. H we have plotted the radial distribution func- 
tions (RDF) of the water oxygen atoms. The minor vari- 
ations among the RDF's of the small methanol system, 
the large methanol system, and the pure water system 
is an indication that the local water structure, as mea- 
sured by this RDF, is at only marginally changed by the 
solvation of a methanol molecule. Note, in this respect, 
that for the 32 molecule the first solvation shell consti- 
tutes a significant fraction of the total number of water 
molecules (see below). 



o.o 
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FIG. 2. Calculated carbon-oxygen (top) and water oxy- 
gen-oxygen (bottom) radial distribution functions for the 
small (solid line) and large (dashed line) methanol system. 

Fig. U also shows the RDF of the methanol carbon and 
water oxygens for the small and large methanol system. 
A pronounced first peak clearly indicates the existence of 



shell of water molecules at a distance of «3.7 A. Com- 
paring the RDF's of the small and large system shows 
a noticeable difference. This should be attributed to 
the limited size of the small system. It suggests that 
a proper description of the solvation structure of a single 
methanol in a cubic periodic simulation box requires at 
least 50 water molecules. Integrating the RDF for the 
large system up to the minimum at r = 5.0 A yields 16 
water molecules in the first solvation shell. The definite 
solvation shell observed in our simulations is consistent 
with the neutron diffraction data of Soper and Finney || 
who studied a 1:9 molar methanol- water system. Differ- 
ences in molarity limits a quantitative comparison of the 
carbon-oxygen RDF, but a qualitative comparison learns 
that peak positions match with the peak values slightly 
more pronounced in the simulation results. 

To analyze the orientational ordering of the water 
molecules around the methanol we computed the distri- 
bution function of the angle between the C-Oh2 bond 
vector and the normal to the plane of the water molecules 
in the first solvation shell. The results show that angle 
distribution is relatively uniform with a small tendency 
towards the tangential orientation, a feature occurs for all 
solvation shell radii in the range of 3.7-5.0 A. Over the 
range of 0°-90° the distribution gradually decays, with 
the value at the tangential orientation (0°) about a factor 
of 2 larger than at the perpendicular orientation (90°). 
Qualitatively, this seems consistent with data for the ori- 
entational distribution obtained from neutron-diffraction 
data |y . However from this experimental data it is con- 
cluded that the water molecules prefer to lie tangential 
and form a cage around the methanol. Our data do not 
give clear evidence for a cage-like structure. However, 
this might be a different interpretation from similar data. 
Note, in this respect, also that the experimental data can- 
not be quantitatively compared to our data, as different 
orientational distribution functions are employed. 

To analyze the hydrogen bonding we adopted the def- 
inition of Ref. 0| : two molecules are hydrogen bonded 
if simultaneously the inter-oxygen distance is less than 
3.5 Aand the OHO angle is smaller than 30°. From 
the simulation of the large system we found that the 
methanol hydroxyl group donates and accepts on average 
0.9 and 1.5 hydrogen bonds, respectively. For a water 
molecule these numbers are equal and measured to be 
1.7 in the simulation of the pure water sample. These 
results indicate that the methanol hydroxyl group par- 
ticipates strongly in the hydrogen bonding network with 
the a donating behavior similar to water hydrogen and a 
accepting character somewhat smaller than a water oxy- 
gen. 

Dynamics 

The time scale (7-10 ps) of the present simulations al- 
lows for a reliable analysis of dynamical properties oc- 
curring on the picosecond time scale. 
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The velocity auto correlation function (VACF) of the 
hydrogen atoms provides an important measure of hy- 
drogen bonding. Fig. |§| shows the Fourier spectrum 
of the calculated VACF of hydrogen atoms of the wa- 
ter molecules in the small and large methanol sample. 
The three distinct peaks correspond to the vibrational 
(3100 cm -1 ), bending (1600 cm -1 ), and librational- 
translational (500 cm -1 ) modes of the water molecules. 
The most important observation is that mutual compar- 
ison of the two methanol samples and the comparison of 
these with the spectrum of the pure water sample (also 
plotted) shows no significant difference, not even for the 
small methanol sample where the solvation shell consti- 
tutes half of the water molecules in the system. This 
demonstrates that also the short-time dynamics of the 
water molecules is hardly affected by the solvation of a 
methanol molecule. 



3.0.2 
0.0 



0.0 





I isolated methanol 
! methanol in 63 water — ■ 




1 31 water + 1 methanol - 

J 63 water + 1 methanol - 

^^^l^^^^^^^j^ 32 water ■ 



500 1000 1500 2000 2500 3000 3500 4000 



V (cm 1 ) 

FIG. 3. Bottom: Calculated Fourier spectrum of the ve- 
locity auto correlation function of the water hydrogens for the 
small methanol system (solid line) , the large methanol system 
(dashed line), and the pure water sample (dotted line). Top: 
Calculated Fourier spectrum of the velocity auto correlation 
function of the hydrogen atom of the methanol hydroxyl group 
for the large methanol system (dashed line) and for an isolated 
methanol molecule (solid line). 

An indication for the average residence time of a water 
molecule in the first solvation shell is obtained by moni- 
toring the trajectories of the individual water molecules. 
We found that in the large methanol system over 7 ps 
10 water molecules left the region within 5 A from the 
methanol carbon. From this we estimate the average res- 
idence time to be of the order of a few picoseconds. 

Fig. U shows the Fourier spectrum of the VACF of the 
hydroxyl H of methanol obtained from the trajectory of 
the large system. The spectrum is of limited accuracy 
due to the relatively short trajectories (7 ps). For com- 
parison, the calculated spectrum for a single methanol 
molecule at T — 200 K is also plotted. In solution the 
OH stretch (^oh) peak, with a calculated gas-phase po- 
sition of about 3300 cm -1 , has shifted by « 200 cm -1 to 
lower frequencies and has a relatively large width. The 
shift and width are both typical characteristics of a hy- 
drogen bond and are also observed in the water spec- 



trum (Fig. H). In contrast to the OH stretch mode, we 
see that the OH-bending mode (<5oh at 1300 cm -1 ) is 
blue-shifted by an amount of 50-100 -1 . A comparison 
with experimental frequency shifts in infrared spectra is 
limited as, to our knowledge, no experimental data for 
dilute methanol-water solutions are reported. However, 
a comparison with measured shifts in liquid methanol 
shows similar trends for the shift of infrared stretch 
(-354 cm -1 ) and bend (+78 cm -1 ) peaks. The torsional 
mode (toh), expected to be shifted upward to around 
600 cm -1 , is not visible in our calculated spectra due to 
the large statistical errors. 

DISCUSSION 

We have studied the solvation of a single methanol 
molecule in water using DFT-based Car-Parrinello 
molecular dynamics simulation. Validation of the ap- 
proach showed that energetics, structural, and dynam- 
ical properties of reference gas-phase compounds were 
sufficient to expect a quantitative accuracy of calculated 
properties. 

The calculated solvation structure supports the exper- 
imental observation ||] that a shell of about 15 water 
molecules is formed around the methanol. Structural 
analysis also learns that the hydrogen bonded network 
of water is only minimally distorted by the presence of 
the methanol molecule. This confirms the proposition 
of Soper et al. Q that speculations that the normal wa- 
ter structure is significantly enhanced by the hydropho- 
bic alkyl group is groundless. The calculations showed 
that methanol OH group is strongly involved in hydro- 
gen bonding, both as acceptor and as donor. Analysis 
of the dynamics learns that the average residence time 
of a water molecule in the first solvation shell is of the 
order of a few picoseconds. The vibrational spectrum of 
the water molecules is hardly changed by the presence 
of the methanol, indicating that the short-time dynam- 
ics is hardly affected by the presence of the methanol 
molecule. Vibrational analysis shows that methanol OH- 
stretch peak is a broad feature that is significantly red- 
shifted upon solvation, confirming its hydrogen-bonding 
character. 

In conclusion, from comparison with available exper- 
imental data we have shown that first-principle DFT- 
based molecular dynamics simulation provides a reason- 
able accurate description of the structure and dynamics 
of a dilute aqueous methanol solution. This opens the 
way towards the study of chemistry involving methanol 
and larger alcohols in water. 
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